The effect of a dominant shrub on annuals plants in a desert system changes along a gradient of precipitation.
Alex Filazzola
May 2017
Ephedra regional gradient (ERG) analyses
Abstract
An under-examined component of the shrub-annual relationship is how regional drivers, such as climate, may alter the sign or magnitude of positive interactions. Interspecific interactions between plants have been shown to be strongly linked to climate, particularly temperature and precipitation. The stress-gradient hypothesis (SGH) predicts that higher abiotic stress (i.e. for deserts, increasing temperature and reduced precipitation) will increase the frequency of positive interactions among shrubs and their annual understory. Regional climate gradients also have indirect effects on plant composition, such as determining consumer abundance and soil nutrient composition. Nutrient availability is particularly affected by precipitation because of altered decomposition rates of organic matter and mineralization. Therefore, the strength of facilitation and operating mechanism of a shrub on the annual plant community may change along a regional gradient.
Hypothesis
We tested the hypothesis that positive interactions among shrubs and annual plants will increase with abiotic stress and reduce nutrient availability along a regional gradient of aridity.
Methods
Weather for growing season
Climate patterns within study
season1.sjd <- season1 %>% filter(Gradient<4) %>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season1.mnp <- season1 %>% filter(Gradient>3)%>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season2.sjd <- season2 %>% filter(Gradient<4) %>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season2.mnp <- season2 %>% filter(Gradient>3)%>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
## Rain vs Temperature in 2016
par(mfrow=c(2,1))
par(mar=c(3.5,4.5,1,4.5))
plot1 <- barplot(height=season1.sjd$Precip, ylim=c(0,14), ylab="Average precipitation San Joaquin (cm)")
points(plot1[,1], season1.sjd$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature at all sites (°C)", 4, line=3)
par(mar=c(4.5,4.5,0,4.5))
plot1 <- barplot(height=season1.mnp$Precip, ylim=c(0,14), ylab="Average precipitation Mojave (cm)")
axis(1, plot1[c(1,30,60,90,120,150,180)], c("Nov","Dec","Jan","Feb","Mar","Apr","May"))
points(plot1[,1], season1.mnp$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature at all sites (°C)", 4, line=3)
## Rain vs Temperature in 2017
par(mfrow=c(2,1))
par(mar=c(3.5,4.5,1,4.5))
plot2 <- barplot(height=season2.sjd$Precip, ylim=c(0,14), ylab="Average precipitation San Joaquin (cm)")
points(plot2[,1], season2.sjd$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature San Joaquin (°C)", 4, line=3)
par(mar=c(4.5,4.5,0,4.5))
plot2 <- barplot(height=season2.mnp$Precip, ylim=c(0,14), ylab="Average precipitation Mojave (cm)")
points(plot2[,1], season2.mnp$min.temp, type="l", col="#FF000050", lwd=2)
axis(1, plot1[c(1,30,60,90,120,150,180)], c("Nov","Dec","Jan","Feb","Mar","Apr","May"))
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature Mojave (°C)", 4, line=3)
### 2016 The rain was inconsistent and mostly absent in the Mojave. This resulted in low germination and producitivty at the southern sites
### 2017 The rain was more plentiful, but in the northern sites, there appears to be a frost period after the majority of the rainfall. Need to check number of frost days
season1.frost <- season1 %>% group_by(Site) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100)
data.frame(season1.frost)## Site frost.days
## 1 Barstow 24.725275
## 2 Cuyama 29.670330
## 3 HeartofMojave 25.824176
## 4 PanocheHills 10.439560
## 5 SheepholeValley 6.043956
## 6 Tecopa 23.626374
## 7 TejonRanch 50.549451season2.frost <- season2 %>% group_by(Site) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100)
data.frame(season2.frost)## Site frost.days
## 1 Barstow 17.679558
## 2 Cuyama 19.889503
## 3 HeartofMojave 16.574586
## 4 PanocheHills 14.917127
## 5 SheepholeValley 1.785714
## 6 Tecopa 25.966851
## 7 TejonRanch 35.911602## Both years had comparable number of frost days
## Compare number of consecutive frost days (i.e. frost periods)
season1[,"frost"] <- ifelse(season1$min.temp<0, -99,season1$min.temp) ## identified days below freezing
season2[,"frost"] <- ifelse(season2$min.temp<0, -99,season2$min.temp) ## identified days below freezing
count.consec <- function(x) {max(rle(as.character(x))$lengths)}
season1.frost <- season1 %>% group_by(Site) %>% summarize(count.consec(frost))
data.frame(season1.frost)## Site count.consec.frost.
## 1 Barstow 10
## 2 Cuyama 10
## 3 HeartofMojave 12
## 4 PanocheHills 5
## 5 SheepholeValley 7
## 6 Tecopa 11
## 7 TejonRanch 14season2.frost <- season2 %>% group_by(Site) %>% summarize(count.consec(frost))
data.frame(season2.frost)## Site count.consec.frost.
## 1 Barstow 5
## 2 Cuyama 6
## 3 HeartofMojave 7
## 4 PanocheHills 6
## 5 SheepholeValley 3
## 6 Tecopa 7
## 7 TejonRanch 13## compare only after plants have germinated
season1.frost <- season1 %>% group_by(Site) %>% filter(year>2015) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100, avg.min.temp=mean(min.temp, na.rm=T))
data.frame(season1.frost)## Site frost.days avg.min.temp
## 1 Barstow 12.396694 5.750413
## 2 Cuyama 11.570248 3.352893
## 3 HeartofMojave 16.528926 4.542149
## 4 PanocheHills 2.479339 6.777686
## 5 SheepholeValley 2.479339 9.394167
## 6 Tecopa 11.570248 6.342149
## 7 TejonRanch 33.884298 1.438017season2.frost <- season2 %>% group_by(Site) %>% filter(year>2016) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100,avg.min.temp=mean(min.temp, na.rm=T))
data.frame(season2.frost)## Site frost.days avg.min.temp
## 1 Barstow 11.666667 6.052500
## 2 Cuyama 16.666667 3.624167
## 3 HeartofMojave 8.333333 5.637838
## 4 PanocheHills 8.333333 6.065833
## 5 SheepholeValley 0.000000 9.386916
## 6 Tecopa 20.833333 5.938333
## 7 TejonRanch 28.333333 2.535833season1.frost <- season1 %>% group_by(Site) %>% filter(year>2015) %>% summarize(count.consec(frost))
data.frame(season1.frost)## Site count.consec.frost.
## 1 Barstow 5
## 2 Cuyama 6
## 3 HeartofMojave 4
## 4 PanocheHills 2
## 5 SheepholeValley 2
## 6 Tecopa 4
## 7 TejonRanch 10season2.frost <- season2 %>% group_by(Site) %>% filter(year>2016)%>% summarize(count.consec(frost))
data.frame(season2.frost)## Site count.consec.frost.
## 1 Barstow 5
## 2 Cuyama 6
## 3 HeartofMojave 3
## 4 PanocheHills 3
## 5 SheepholeValley 2
## 6 Tecopa 5
## 7 TejonRanch 10Microenvironmental differences
We compared temperature and relative humidity between shrub and open microsites among all sites along the regional gradient 
Shapiro-Wilk normality test
data: fit1$residuals
W = 0.87881, p-value < 2.2e-16 Df Sum Sq Mean Sq F value Pr(>F)
Microsite 1 26.5 26.51 59.886 1.24e-14 ***
gradient 6 233.9 38.98 88.065 < 2e-16 ***
Microsite:gradient 6 13.1 2.18 4.927 4.95e-05 ***
Residuals 4290 1898.9 0.44
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingnessAnalysis of Deviance Table
Model: binomial, link: logit
Response: RH/100
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 4333 974.99
Microsite 1 1.33 4332 973.66 0.2492
gradient 6 528.18 4326 445.48 <2e-16 ***
Microsite:gradient 6 7.69 4320 437.79 0.2613
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Analysis of Deviance Table
Model: Gamma, link: inverse
Response: swc
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev F Pr(>F)
NULL 404 231.405
Site 6 189.683 398 41.722 381.2179 <2e-16 ***
Microsite 1 0.287 397 41.435 3.4579 0.0637 .
Site:Microsite 6 0.696 391 40.739 1.3994 0.2136
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Analysis of Deviance Table
Model: Gamma, link: inverse
Response: swc
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev F Pr(>F)
NULL 419 183.193
Site 6 144.930 413 38.263 279.4537 <2e-16 ***
Microsite 1 0.080 412 38.183 0.9284 0.3358
Site:Microsite 6 0.535 406 37.647 1.0322 0.4037
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
site level means
## site coordintes
site.gps <- read.csv("Data/ERGsites.csv")
# nutrient data for each site
nutrients <- read.csv("Data/ERG.soilnutrients.csv")
## long term climate data extracted from worldclim
clim.dat <- read.csv("Data/ERG.worldclim.csv")
## check correlation among long-term climate variables
cor(clim.dat[,3:8]) ## Temp.wettest.QR least correlated ## Annual.Temp Temp.Seasonality Temp.wettest.QR
## Annual.Temp 1.0000000 0.9051791 0.7260617
## Temp.Seasonality 0.9051791 1.0000000 0.4650548
## Temp.wettest.QR 0.7260617 0.4650548 1.0000000
## Annual.precip -0.8911069 -0.9650533 -0.5313460
## Precip.seasonality -0.9715419 -0.9774646 -0.6165170
## Precipt.wettest.QR -0.8839124 -0.9570957 -0.5523347
## Annual.precip Precip.seasonality Precipt.wettest.QR
## Annual.Temp -0.8911069 -0.9715419 -0.8839124
## Temp.Seasonality -0.9650533 -0.9774646 -0.9570957
## Temp.wettest.QR -0.5313460 -0.6165170 -0.5523347
## Annual.precip 1.0000000 0.9613729 0.9980959
## Precip.seasonality 0.9613729 1.0000000 0.9543657
## Precipt.wettest.QR 0.9980959 0.9543657 1.0000000## Obtain mean shrub traits for each site
shrubs <- read.csv("Data/ERG.shrub.csv")
shrubs <- subset(shrubs, Microsite=="shrub")
shrubs.mean <- shrubs %>% group_by(Gradient,Site) %>% summarise_all(funs(mean))
shrubs.mean <- data.frame(shrubs.mean)
shrubs.vars <- shrubs.mean[,c("Site","Gradient","volume","canopy","Dx","DxEph","Compaction")]
## test correlation among shrub traits
cor(shrubs.vars[,3:7]) ## compaction x volume & DxEph x Dx high correlation## volume canopy Dx DxEph Compaction
## volume 1.0000000 0.2243902 0.40036651 0.2631563 -0.64927168
## canopy 0.2243902 1.0000000 -0.68479570 -0.3538094 -0.20833287
## Dx 0.4003665 -0.6847957 1.00000000 0.7803527 -0.06531485
## DxEph 0.2631563 -0.3538094 0.78035272 1.0000000 0.33769574
## Compaction -0.6492717 -0.2083329 -0.06531485 0.3376957 1.00000000vifstep(shrubs.vars[,3:7], th=10) ## Dx showing collinearity problems## 1 variables from the 5 input variables have collinearity problem:
##
## Dx
##
## After excluding the collinear variables, the linear correlation coefficients ranges between:
## min correlation ( Compaction ~ canopy ): -0.2083329
## max correlation ( Compaction ~ volume ): -0.6492717
##
## ---------- VIFs of the remained variables --------
## Variables VIF
## 1 volume 4.043473
## 2 canopy 1.474526
## 3 DxEph 2.874621
## 4 Compaction 3.760284shrub.site <- shrubs.vars[,-c(1,2,5)]
rownames(shrub.site) <- shrubs.vars[,1]
## PCA of shrub characteristics
pca1 <- prcomp(log(shrub.site), scale=T)
plot(pca1)
biplot(pca1, scale = T)
summary(pca1) ## 77% variation explained## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 1.2958 1.1879 0.8922 0.33722
## Proportion of Variance 0.4198 0.3528 0.1990 0.02843
## Cumulative Proportion 0.4198 0.7725 0.9716 1.00000## PCA of site characteristics for season1
## Obtain mean weather variables for each site
season1.mean <- season1 %>% group_by(Gradient,Site) %>% summarise(temp.var=var(avg.temp,na.rm=T),Precip=sum(Precip),wind=mean(wind.speed, na.rm=T))
season1.mean <- data.frame(season1.mean)
## extract key variables
## dropped RH, and chose min.temp because least correlation with others & cold stress
##combine nutrients and long-term averages
site.vars <- data.frame(season1.mean[,3:5],site.gps["elevation"]) ## drop Phosphorus and other climate variables because of correlations,
row.names(site.vars) <- shrubs.vars[,1]
cor(site.vars)## temp.var Precip wind elevation
## temp.var 1.00000000 -0.8679054 0.05816084 -0.2699588
## Precip -0.86790537 1.0000000 -0.37279135 0.1222504
## wind 0.05816084 -0.3727914 1.00000000 0.1054497
## elevation -0.26995879 0.1222504 0.10544974 1.0000000site.vars2016 <- site.vars
## check for collinearity
vifstep(site.vars, th=10) ## remove potassium and temperature minimum## No variable from the 4 input variables has collinearity problem.
##
## The linear correlation coefficients ranges between:
## min correlation ( wind ~ temp.var ): 0.05816084
## max correlation ( Precip ~ temp.var ): -0.8679054
##
## ---------- VIFs of the remained variables --------
## Variables VIF
## 1 temp.var 6.640710
## 2 Precip 7.316232
## 3 wind 1.739558
## 4 elevation 1.142681pca1 <- prcomp(log(abs(site.vars)), scale=T)
plot(pca1)
biplot(pca1)
## check contribution of loadings
pca1$rotation## PC1 PC2 PC3 PC4
## temp.var 0.6001002 -0.1669018 0.50554765 -0.5970303
## Precip -0.6703313 -0.1105034 -0.09936844 -0.7270288
## wind 0.2743589 0.7981791 -0.43410209 -0.3149488
## elevation -0.3395039 0.5681927 0.73898773 0.1256633aload <- abs(pca1$rotation)
sweep(aload, 2, colSums(aload), "/")## PC1 PC2 PC3 PC4
## temp.var 0.3184748 0.1015355 0.28433406 0.33832382
## Precip 0.3557466 0.0672253 0.05588758 0.41199111
## wind 0.1456030 0.4855763 0.24415109 0.17847449
## elevation 0.1801756 0.3456629 0.41562726 0.07121058summary(pca1) ## 89% variation explained## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 1.456 1.0458 0.8577 0.22479
## Proportion of Variance 0.530 0.2734 0.1839 0.01263
## Cumulative Proportion 0.530 0.8034 0.9874 1.00000gradient1.season1 <- pca1$x[,1]
gradient2.season1 <- pca1$x[,2]
## season2
season2.mean <- season2 %>% group_by(Gradient,Site) %>% summarise(temp.var=var(avg.temp,na.rm=T),Precip=sum(Precip, na.rm=T),wind=mean(wind.speed, na.rm=T))
season2.mean <- data.frame(season2.mean)
## extract key variables
## dropped RH, and chose min.temp because least correlation with others & cold stress
##combine nutrients and long-term averages
site.vars <- data.frame(season2.mean[,3:5],site.gps["elevation"]) ## drop Phosphorus and other climate variables because of correlations,
row.names(site.vars) <- shrubs.vars[,1]
cor(site.vars)## temp.var Precip wind elevation
## temp.var 1.0000000 -0.7427028 -0.2262409 -0.3715471
## Precip -0.7427028 1.0000000 -0.1551517 0.1813600
## wind -0.2262409 -0.1551517 1.0000000 0.1052805
## elevation -0.3715471 0.1813600 0.1052805 1.0000000site.vars2017 <- site.vars
## check for collinearity
vifstep(site.vars, th=10) ## remove potassium and temperature minimum## No variable from the 4 input variables has collinearity problem.
##
## The linear correlation coefficients ranges between:
## min correlation ( elevation ~ wind ): 0.1052805
## max correlation ( Precip ~ temp.var ): -0.7427028
##
## ---------- VIFs of the remained variables --------
## Variables VIF
## 1 temp.var 3.438786
## 2 Precip 3.035337
## 3 wind 1.402010
## 4 elevation 1.192010pca2 <- prcomp(log(abs(site.vars)), scale=T)
plot(pca2)
biplot(pca2)
summary(pca2) ## 85% variation explained## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 1.434 1.0121 0.8809 0.37883
## Proportion of Variance 0.514 0.2561 0.1940 0.03588
## Cumulative Proportion 0.514 0.7701 0.9641 1.00000## check contribution of loadings
pca2$rotation## PC1 PC2 PC3 PC4
## temp.var 0.6557532 -0.09123357 0.20565681 0.7206729
## Precip -0.6176537 -0.16891771 -0.40105332 0.6550778
## wind -0.1019884 0.97279862 0.06826009 0.1964733
## elevation -0.4220070 -0.12963831 0.89005734 0.1135866aload <- abs(pca2$rotation)
sweep(aload, 2, colSums(aload), "/")## PC1 PC2 PC3 PC4
## temp.var 0.36483385 0.06695609 0.1314078 0.42749339
## Precip 0.34363685 0.12396828 0.2562596 0.38858328
## wind 0.05674213 0.71393442 0.0436159 0.11654531
## elevation 0.23478717 0.09514122 0.5687167 0.06737801## define gradients
gradient1.season2 <- pca2$x[,1]
gradient2.season2 <- pca2$x[,2]Compare gradient to phytometer
# crs.world <-CRS("+proj=longlat +datum=WGS84")
# gps <- data.frame(x=site.gps$long,y=site.gps$lat)
# coordinates(gps) <- ~x+y
# proj4string(gps) <- crs.world
#
# ## Download and extract PET/aridity
# r1 <- raster("C:\\Users\\Alessandro\\Downloads\\Global Aridity - Annual\\AI_annual\\ai_yr\\w001001.adf", package="raster") #aridity
# r2 <- raster("C:\\Users\\Alessandro\\Downloads\\Global PET - Annual\\PET_he_annual\\pet_he_yr\\w001001.adf", package="raster") #potential evapotranspiration
# calclim <- stack(r1,r2)
# names(calclim) <- c("aridity","PET")
# arid.vals <- raster::extract(calclim, gps)
# rownames(arid.vals) <- site.gps[,"name"]
# colnames(arid.vals)[1] <- "Gradient"
# write.csv(arid.vals, "Data/aridity.PET.csv")
arid.vals <- read.csv("Data/aridity.PET.csv")
colnames(arid.vals)[1] <- "Site"
par(mfrow=c(1,2))
mean.phyto <- census %>% filter(Census=="end") %>%group_by(Year, Gradient, Site, Microsite) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
mean.phyto <- data.frame(mean.phyto)
rii.data <- rii(census, 1:6, var=9:15)
rii.mean <- rii.data %>% group_by(Year, Gradient,Site) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
rii.mean <- data.frame(rii.mean)
season1.mean <- season1 %>% group_by(Gradient,Site) %>% summarise(temp.var=var(avg.temp,na.rm=T),Precip=sum(Precip),wind=mean(wind.speed, na.rm=T), max.temp=abs(mean(max.temp, na.rm=T)),min.temp=abs(mean(min.temp, na.rm=T)),avg.temp=abs(mean(avg.temp, na.rm=T)))
s1.mean <- data.frame(season1.mean)
season2.mean <- season2 %>% group_by(Gradient,Site) %>% summarise(temp.var=var(avg.temp,na.rm=T),Precip=sum(Precip, na.rm=T),wind=mean(wind.speed, na.rm=T), max.temp=abs(mean(max.temp, na.rm=T)),min.temp=abs(mean(min.temp, na.rm=T)),avg.temp=abs(mean(avg.temp, na.rm=T)))
s2.mean <- data.frame(season2.mean)
s1.mean[,"arid.gradient"] <- log(s1.mean[,"Precip"]/arid.vals[,"PET"])
s2.mean[,"arid.gradient"] <- log(s2.mean[,"Precip"]/(arid.vals[,"PET"]))
## collect aridity gradient data
site.climate <- rbind(s1.mean,s2.mean)
site.climate[,"Year"] <- as.factor(c(rep("2016",7),rep("2017",7)))
## combine climate data with census
mean.phyto <- merge(mean.phyto , site.climate, by=c("Year","Site"))
phyto.2016 <- subset(mean.phyto, Year==2016)
plot(jitter(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"],1), phyto.2016[phyto.2016$Microsite=="open","phyto.biomass"], ylim=c(0,2), ylab=c("all phytometer biomass (g)"), xlab=c("aridity gradient"))
points(jitter(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"],2), phyto.2016[phyto.2016$Microsite=="shrub","phyto.biomass"], pch=19)
legend(-4.5,1.8, c("shrub","open"), pch=c(19,1))
text(-2.3, 1.8, "2016", cex=1.5)
plot(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"], rii.mean[rii.mean$Year==2016,"phyto.biomass"], ylim=c(-0.4,0.4), ylab="RII all phytometer biomass", pch=19, xlab=c("aridity gradient"))
text(-2.3, 0.38, "2016", cex=1.5)
phyto.2017 <- subset(mean.phyto, Year==2017)
plot(jitter(phyto.2017[phyto.2017$Microsite=="open","arid.gradient"],1), phyto.2017[phyto.2017$Microsite=="open","phyto.biomass"], ylim=c(0,3), xlab=c("aridity gradient"),ylab=c("all phytometer biomass (g)"))
points(jitter(phyto.2017[phyto.2017$Microsite=="open","arid.gradient"],2), phyto.2017[phyto.2017$Microsite=="shrub","phyto.biomass"], pch=19)
legend(-4.5,1.8, c("shrub","open"), pch=c(19,1))
text(-2.0, 2.8, "2017", cex=1.5)
plot(phyto.2017[phyto.2017$Microsite=="open","arid.gradient"], rii.mean[rii.mean$Year==2017,"phyto.biomass"], ylim=c(-0.4,0.4), ylab="RII all phytometer biomass", pch=19, xlab=c("aridity gradient"))
text(-2.0, 0.38, "2017", cex=1.5)
## regress gradient and microsite on biomass for phytometers in 2016
m1 <- lm(ihs(Phacelia.biomass)~arid.gradient * Microsite, data=phyto.2016)
summary(m1) ## nothing significant##
## Call:
## lm(formula = ihs(Phacelia.biomass) ~ arid.gradient * Microsite,
## data = phyto.2016)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.26646 -0.16692 0.01821 0.09053 0.43991
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.75896 0.31899 2.379 0.0387 *
## arid.gradient 0.16726 0.09147 1.829 0.0974 .
## Micrositeshrub 0.81930 0.45112 1.816 0.0994 .
## arid.gradient:Micrositeshrub 0.19010 0.12936 1.470 0.1724
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2383 on 10 degrees of freedom
## Multiple R-squared: 0.674, Adjusted R-squared: 0.5762
## F-statistic: 6.893 on 3 and 10 DF, p-value: 0.008498m2 <- lm(ihs(Plantago.biomass)~arid.gradient * Microsite, data=phyto.2016)
summary(m2) ## nothing significant##
## Call:
## lm(formula = ihs(Plantago.biomass) ~ arid.gradient * Microsite,
## data = phyto.2016)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.070995 -0.048467 0.003316 0.024049 0.115472
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.27085 0.07901 3.428 0.00646 **
## arid.gradient 0.05886 0.02266 2.598 0.02658 *
## Micrositeshrub -0.12531 0.11174 -1.121 0.28831
## arid.gradient:Micrositeshrub -0.02928 0.03204 -0.914 0.38239
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05904 on 10 degrees of freedom
## Multiple R-squared: 0.4794, Adjusted R-squared: 0.3232
## F-statistic: 3.069 on 3 and 10 DF, p-value: 0.07774m3 <- lm(ihs(Salvia.biomass)~arid.gradient * Microsite, data=phyto.2016)
summary(m3)## aridity significant##
## Call:
## lm(formula = ihs(Salvia.biomass) ~ arid.gradient * Microsite,
## data = phyto.2016)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.15051 -0.09099 -0.01568 0.03332 0.26045
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.748316 0.194293 3.851 0.0032 **
## arid.gradient 0.159284 0.055714 2.859 0.0170 *
## Micrositeshrub -0.000770 0.274771 -0.003 0.9978
## arid.gradient:Micrositeshrub 0.004763 0.078792 0.060 0.9530
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1452 on 10 degrees of freedom
## Multiple R-squared: 0.6281, Adjusted R-squared: 0.5165
## F-statistic: 5.63 on 3 and 10 DF, p-value: 0.01599## regress gradient and microsite on biomass for phytometers in 2017
m4 <- lm(ihs(Phacelia.biomass)~arid.gradient * Microsite, data=phyto.2017)
summary(m4) ## nothing significant##
## Call:
## lm(formula = ihs(Phacelia.biomass) ~ arid.gradient * Microsite,
## data = phyto.2017)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.53579 -0.12720 0.00415 0.29889 0.40793
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.725 3.895 -1.727 0.135
## arid.gradient -2.514 1.315 -1.912 0.104
## Micrositeshrub 6.929 4.012 1.727 0.135
## arid.gradient:Micrositeshrub 2.341 1.363 1.718 0.137
##
## Residual standard error: 0.4168 on 6 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.3953, Adjusted R-squared: 0.09289
## F-statistic: 1.307 on 3 and 6 DF, p-value: 0.3557m5 <- lm(ihs(Plantago.biomass)~arid.gradient * Microsite, data=phyto.2017)
summary(m5) ## aridity significant##
## Call:
## lm(formula = ihs(Plantago.biomass) ~ arid.gradient * Microsite,
## data = phyto.2017)
##
## Residuals:
## 15 16 19 20 23 24 25
## -0.008390 -0.070144 0.006678 -0.031325 0.126617 0.249998 -0.124905
## 26
## -0.148529
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.4847 1.6331 -0.909 0.415
## arid.gradient -0.6486 0.5513 -1.177 0.305
## Micrositeshrub 3.1150 2.3096 1.349 0.249
## arid.gradient:Micrositeshrub 1.1087 0.7796 1.422 0.228
##
## Residual standard error: 0.1748 on 4 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.4909, Adjusted R-squared: 0.109
## F-statistic: 1.285 on 3 and 4 DF, p-value: 0.3935m6 <- lm(ihs(Salvia.biomass)~arid.gradient * Microsite, data=phyto.2017)
summary(m6)## nothing significant##
## Call:
## lm(formula = ihs(Salvia.biomass) ~ arid.gradient * Microsite,
## data = phyto.2017)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.52276 -0.13176 -0.02786 0.05787 0.53706
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5442 0.6266 -0.869 0.405
## arid.gradient -0.4200 0.2423 -1.733 0.114
## Micrositeshrub 0.2852 0.8861 0.322 0.754
## arid.gradient:Micrositeshrub 0.1488 0.3427 0.434 0.673
##
## Residual standard error: 0.335 on 10 degrees of freedom
## Multiple R-squared: 0.3112, Adjusted R-squared: 0.1046
## F-statistic: 1.506 on 3 and 10 DF, p-value: 0.2723#
# ## Phacelia season 1
# plot(jitter(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"] ,1), phyto.2016[phyto.2016$Microsite=="open","Phacelia.biomass"], ylim=c(0,2))
# points(jitter(phyto.2016[phyto.2016$Microsite=="shrub","arid.gradient"] ,2), phyto.2016[phyto.2016$Microsite=="shrub","Phacelia.biomass"], pch=19)
#
# m1 <- lm(ihs(Phacelia.biomass)~arid.gradient * Microsite, data=phyto.2016)
# m2 <- lm(ihs(Plantago.biomass)~arid.gradient * Microsite, data=phyto.2016)
# m3 <- lm(ihs(Salvia.biomass)~arid.gradient * Microsite, data=phyto.2016)
#
# plot(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"] , rii.mean[rii.mean$Year==2016,"Phacelia.biomass"], ylim=c(-0.4,0.4))
#
# ## Plantago season 1
# plot(jitter(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"] ,1), phyto.2016[phyto.2016$Microsite=="open","Plantago.biomass"], ylim=c(0,2))
# points(jitter(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"] ,2), phyto.2016[phyto.2016$Microsite=="shrub","Plantago.biomass"], pch=19)
#
# plot(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"] , rii.mean[rii.mean$Year==2016,"Plantago.biomass"], ylim=c(-0.4,0.4))
#
# ## Salvia season 1
# plot(jitter(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"] ,1), phyto.2016[phyto.2016$Microsite=="open","Salvia.biomass"], ylim=c(0,2))
# points(jitter(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"] ,2), phyto.2016[phyto.2016$Microsite=="shrub","Salvia.biomass"], pch=19)
#
# plot(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"] , rii.mean[rii.mean$Year==2016,"Salvia.biomass"], ylim=c(-0.4,0.4))
#
#
#
# ## Phacelia season 2
# plot(jitter(gradient1.season2,1), phyto.2017[phyto.2017$Microsite=="open","Phacelia.biomass"], ylim=c(0,2))
# points(jitter(gradient1.season2,2), phyto.2017[phyto.2017$Microsite=="shrub","Phacelia.biomass"], pch=19)
#
# plot(gradient1.season2, rii.mean[rii.mean$Year==2017,"Phacelia.biomass"], ylim=c(-0.4,0.4))
#
#
# ## Plantago season 2
# plot(jitter(gradient1.season2,1), phyto.2017[phyto.2017$Microsite=="open","Plantago.biomass"], ylim=c(0,2))
# points(jitter(gradient1.season2,2), phyto.2017[phyto.2017$Microsite=="shrub","Plantago.biomass"], pch=19)
#
# plot(gradient1.season2, rii.mean[rii.mean$Year==2017,"Plantago.biomass"], ylim=c(-0.4,0.4))
#
# ## Salvia season 2
# plot(jitter(gradient1.season2,1), phyto.2017[phyto.2017$Microsite=="open","Salvia.biomass"], ylim=c(0,2))
# points(jitter(gradient1.season2,2), phyto.2017[phyto.2017$Microsite=="shrub","Salvia.biomass"], pch=19)
#
#
# plot(gradient1.season2, rii.mean[rii.mean$Year==2017,"Salvia.biomass"], ylim=c(-0.4,0.4))
## both seasons
## responses
plot(mean.phyto[mean.phyto$Microsite=="open","arid.gradient"], mean.phyto[mean.phyto$Microsite=="open","phyto.biomass"], ylim=c(0,3), ylab=c("all phytometer biomass (g)"))
points(jitter(mean.phyto[mean.phyto$Microsite=="shrub","arid.gradient"],1), mean.phyto[mean.phyto$Microsite=="shrub","phyto.biomass"], pch=19)
legend(-4.5,2.8, c("shrub","open"), pch=c(19,1))
text(-2.3, 2.8, "both years", cex=1.3)
#
# ## Phacelia
# plot(jitter(mean.phyto[mean.phyto$Microsite=="open","arid.gradient"],1), mean.phyto[mean.phyto$Microsite=="open","Phacelia.biomass"], ylim=c(0,2))
# points(jitter(mean.phyto[mean.phyto$Microsite=="shrub","arid.gradient"],2), mean.phyto[mean.phyto$Microsite=="shrub","Phacelia.biomass"], pch=19)
# m1 <- lm(ihs(Phacelia.biomass) ~ arid.gradient * Microsite, data=mean.phyto)
# summary(m1)
#
# ## Plantago
# plot(jitter(mean.phyto[mean.phyto$Microsite=="open","arid.gradient"],1), mean.phyto[mean.phyto$Microsite=="open","Plantago.biomass"], ylim=c(0,0.5))
# points(jitter(mean.phyto[mean.phyto$Microsite=="shrub","arid.gradient"],2), mean.phyto[mean.phyto$Microsite=="shrub","Plantago.biomass"], pch=19)
# m2 <- lm(ihs(Phacelia.biomass) ~ arid.gradient * Microsite, data=mean.phyto)
# summary(m2)
#
#
# ## Salvia season 1
# plot(jitter(mean.phyto[mean.phyto$Microsite=="open","arid.gradient"],1), mean.phyto[mean.phyto$Microsite=="open","Salvia.biomass"], ylim=c(0,2))
# points(jitter(mean.phyto[mean.phyto$Microsite=="shrub","arid.gradient"],2), mean.phyto[mean.phyto$Microsite=="shrub","Salvia.biomass"], pch=19)
# m3 <- lm(ihs(Salvia.biomass) ~ arid.gradient * Microsite, data=mean.phyto)
# summary(m3)
### Rii plots
rii.mean <- merge(rii.mean, site.climate, by=c("Year","Site"))
## biomass by RII
plot(mean.phyto[mean.phyto$Microsite=="open","Phacelia.biomass"], rii.mean[,"Phacelia.biomass"], ylim=c(-0.3,0.3), pch=19, ylab="RII biomass", xlab="average phytometer biomass per site")
points(mean.phyto[mean.phyto$Microsite=="open","Plantago.biomass"], rii.mean[,"Plantago.biomass"], pch=21, bg="blue")
points(mean.phyto[mean.phyto$Microsite=="open","Salvia.biomass"], rii.mean[,"Salvia.biomass"], pch=21, bg="red")
abline(h=0, lty=2, lwd=2)
legend(0.38,-0.15, c("Phacelia","Plantago","Salvia"), pch=21, pt.bg=c("black","blue","red"))
## rii variability
plot(rii.mean[,"arid.gradient"], rii.mean[,"Phacelia.biomass"], ylim=c(-0.5,0.5), pch=19, ylab="Rii phytometer biomass", xlab="aridity gradient")
points(rii.mean[,"arid.gradient"], rii.mean[,"Plantago.biomass"], pch=21, bg="blue")
points(rii.mean[,"arid.gradient"], rii.mean[,"Salvia.biomass"], pch=21, bg="red")
abline(h=0, lty=2, lwd=2)
legend(-4.5, 0.45, c("Phacelia","Plantago","Salvia"), pch=21, pt.bg=c("black","blue","red"))
## compare convex hull
## hull function
Plot_ConvexHull<-function(xcoord, ycoord, lcolor){
hpts <- chull(x = xcoord, y = ycoord)
hpts <- c(hpts, hpts[1])
lines(xcoord[hpts], ycoord[hpts], col = lcolor)
}
Plot_ConvexHull(xcoord = rii.mean[,"arid.gradient"], ycoord = rii.mean[,"Phacelia.biomass"], lcolor = "black")
Plot_ConvexHull(xcoord = rii.mean[,"arid.gradient"], ycoord = rii.mean[,"Plantago.biomass"], lcolor = "blue")
Plot_ConvexHull(xcoord = rii.mean[,"arid.gradient"], ycoord = rii.mean[,"Salvia.biomass"], lcolor = "red")
r.df <- rii.mean[,c("arid.gradient","Phacelia.biomass","Plantago.biomass","Salvia.biomass")]
## transform data into long format
r.df <- gather(r.df, species, biomass, 2:4)
## function to create a convex hull around data
hull.time <- function(df){
df[chull(df$arid.gradient,df$biomass),] # chull really is useful, even outside of contrived examples.
}
## break data into lists and recombined to apply polygon
splitData <- split(r.df, r.df$species)
appliedData <- lapply(splitData, hull.time)
combinedData <- do.call(rbind, appliedData)
## create polygon
poly1 <- ggplot(r.df, aes(x=arid.gradient, y=biomass))+
geom_polygon(data = combinedData, # This is also a nice example of how to plot
aes(x=arid.gradient, y=biomass, group = species, color=species), fill=NA, # two superimposed geoms
alpha = 1/2) + theme_bw()
poly1
## regress rii on gradient
m1 <- lm(Phacelia.biomass~arid.gradient, data=rii.mean)
summary(m1) ## not significant##
## Call:
## lm(formula = Phacelia.biomass ~ arid.gradient, data = rii.mean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.022533 -0.013199 -0.005837 0.008594 0.058681
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.038659 0.020960 1.844 0.0899 .
## arid.gradient 0.008913 0.006828 1.305 0.2162
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02266 on 12 degrees of freedom
## Multiple R-squared: 0.1244, Adjusted R-squared: 0.05139
## F-statistic: 1.704 on 1 and 12 DF, p-value: 0.2162m2 <- lm(Plantago.biomass~arid.gradient, data=rii.mean)
summary(m2) ## not significant##
## Call:
## lm(formula = Plantago.biomass ~ arid.gradient, data = rii.mean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.030560 -0.002438 0.002575 0.006742 0.019216
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.006458 0.014218 -0.454 0.658
## arid.gradient -0.001721 0.004632 -0.372 0.717
##
## Residual standard error: 0.01537 on 12 degrees of freedom
## Multiple R-squared: 0.01137, Adjusted R-squared: -0.07101
## F-statistic: 0.1381 on 1 and 12 DF, p-value: 0.7167m3 <- lm(Salvia.biomass~arid.gradient, data=rii.mean)
summary(m3) ## not significant##
## Call:
## lm(formula = Salvia.biomass ~ arid.gradient, data = rii.mean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.051586 -0.015918 -0.001203 0.007707 0.066895
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.015346 0.029966 -0.512 0.618
## arid.gradient -0.004097 0.009762 -0.420 0.682
##
## Residual standard error: 0.0324 on 12 degrees of freedom
## Multiple R-squared: 0.01447, Adjusted R-squared: -0.06766
## F-statistic: 0.1761 on 1 and 12 DF, p-value: 0.6821#
# m1 <- lm(Phacelia.biomass~arid.gradient, data=subset(rii.mean,Year == 2016))
# m2 <- lm(Plantago.biomass~arid.gradient, data=subset(rii.mean,Year == 2016))
# m3 <- lm(Salvia.biomass~arid.gradient, data=subset(rii.mean,Year == 2016))
#
#
# m1 <- lm(Phacelia.biomass~arid.gradient, data=subset(rii.mean,Year == 2017))
# m2 <- lm(Plantago.biomass~arid.gradient, data=subset(rii.mean,Year == 2017))
# m3 <- lm(Salvia.biomass~arid.gradient, data=subset(rii.mean,Year == 2017))
#
#
# plot(rii.mean[,"min.temp"], rii.mean[,"Phacelia.biomass"], ylim=c(-0.4,0.4), pch=19)
# points(rii.mean[,"min.temp"], rii.mean[,"Plantago.biomass"], pch=21, bg="blue")
# points(rii.mean[,"min.temp"], rii.mean[,"Salvia.biomass"], pch=21, bg="red")
# abline(h=0, lty=2, lwd=2)
#
# rii.test <- subset(rii.mean, Gradient.x>3 & Year == 2017 | Year == 2016)
#
#
# plot(rii.test[,"arid.gradient"], rii.test[,"Phacelia.biomass"], ylim=c(-0.4,0.4), pch=19)
# points(rii.test[,"arid.gradient"], rii.test[,"Plantago.biomass"], pch=21, bg="blue")
# points(rii.test[,"arid.gradient"], rii.test[,"Salvia.biomass"], pch=21, bg="red")
# abline(h=0, lty=2, lwd=2)
#
#
# plot(rii.mean[,"arid.gradient"], rii.mean[,"Phacelia"], ylim=c(-0.4,0.4), pch=19)
# points(rii.mean[,"arid.gradient"], rii.mean[,"Plantago"], pch=21, bg="blue")
# points(rii.mean[,"arid.gradient"], rii.mean[,"Salvia"], pch=21, bg="red")
# abline(h=0, lty=2, lwd=2)
#
#
# m1 <- lm(Phacelia~ poly(arid.gradient,2), data=rii.mean)
# summary(m1)
# m2 <- lm(Plantago.biomass~arid.gradient, data=rii.mean)
# m3 <- lm(Salvia.biomass~arid.gradient, data=rii.mean)
Differences from Nutrient content
## join aridity with nutrients
nutrients.climate <- merge(nutrients,subset(site.climate, Year==2016), by=c("Site"))
## Nitrogen difference between shrub and sites
m.nit <- aov(log(N) ~ Site * microsite, data=nutrients)
summary(m.nit)## Df Sum Sq Mean Sq F value Pr(>F)
## Site 6 63.04 10.51 13.241 2.86e-09 ***
## microsite 1 40.16 40.16 50.614 2.24e-09 ***
## Site:microsite 6 4.96 0.83 1.042 0.408
## Residuals 56 44.43 0.79
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(m.nit, "microsite")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = log(N) ~ Site * microsite, data = nutrients)
##
## $microsite
## diff lwr upr p adj
## shrub-open 1.514873 1.088317 1.941429 0## Potassium difference between shrub and sites
m.pot <- aov(log(K) ~ Site * microsite, data=nutrients)
summary(m.pot)## Df Sum Sq Mean Sq F value Pr(>F)
## Site 6 15.897 2.650 22.013 4.09e-13 ***
## microsite 1 4.445 4.445 36.934 1.14e-07 ***
## Site:microsite 6 1.605 0.268 2.223 0.0541 .
## Residuals 56 6.740 0.120
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(m.pot, "microsite")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = log(K) ~ Site * microsite, data = nutrients)
##
## $microsite
## diff lwr upr p adj
## shrub-open 0.5040081 0.3378755 0.6701406 1e-07## Phosphorus difference between shrub and sites
m.pho <- aov(log(P) ~ Site * microsite, data=nutrients)
summary(m.pho)## Df Sum Sq Mean Sq F value Pr(>F)
## Site 6 23.163 3.861 14.33 8.10e-10 ***
## microsite 1 10.546 10.546 39.15 5.79e-08 ***
## Site:microsite 6 3.394 0.566 2.10 0.0676 .
## Residuals 56 15.085 0.269
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(m.pho, "microsite")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = log(P) ~ Site * microsite, data = nutrients)
##
## $microsite
## diff lwr upr p adj
## shrub-open 0.7762734 0.5277317 1.024815 1e-07nutrient.mean <- nutrients.climate %>% group_by(Site, arid.gradient) %>% summarize(nitrogen=mean(N), phosphorus=mean(P), potassium=mean(K))
nutrient.mean <- data.frame(nutrient.mean)
ggplot(nutrient.mean) + geom_jitter(aes(x=arid.gradient, y=nitrogen), color="#E69F00", size=3) + stat_smooth(method="lm", formula= y~poly(x,2),aes(x=arid.gradient, y=nitrogen), color="#E69F00") + geom_jitter(aes(x=arid.gradient, y=potassium/20), color="#56B4E9", size=3) + stat_smooth(method="lm", formula= y~poly(x,2),aes(x=arid.gradient, y=potassium/20), color="#56B4E9") + geom_jitter(aes(x=arid.gradient, y=phosphorus), color="#009E73", size=3) +ylab("soil nutrient content")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab("aridity gradient")+
annotate("segment", x = -4.5, xend = -4.35, y = 25, yend = 25, colour = "#E69F00", size=2) + ## add custom legend
annotate("text", x = -4.3, y = 25, label = "Nitrogen", size=5, hjust=0)+
annotate("segment", x = -4.5, xend = -4.35, y = 23, yend = 23, colour = "#56B4E9", size=2) + ## add custom legend
annotate("text", x = -4.3, y = 23, label = "Potassium", size=5, hjust=0)+
annotate("segment", x = -4.5, xend = -4.35, y = 21, yend = 21, colour = "#009E73", size=2) + ## add custom legend
annotate("text", x = -4.3, y = 21, label = "Phosphorus", size=5, hjust=0)
m1 <- lm(nitrogen~poly(arid.gradient,2), data=nutrient.mean)
summary(m1)##
## Call:
## lm(formula = nitrogen ~ poly(arid.gradient, 2), data = nutrient.mean)
##
## Residuals:
## 1 2 3 4 5 6 7
## 1.1092 1.0570 -0.4565 -0.8511 -0.8789 -0.3757 0.3960
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.4326 0.3942 13.782 0.000161 ***
## poly(arid.gradient, 2)1 -12.6364 1.0429 -12.116 0.000266 ***
## poly(arid.gradient, 2)2 7.0832 1.0429 6.792 0.002454 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.043 on 4 degrees of freedom
## Multiple R-squared: 0.9797, Adjusted R-squared: 0.9695
## F-statistic: 96.47 on 2 and 4 DF, p-value: 0.0004126m2 <- lm(phosphorus~arid.gradient, data=nutrient.mean)
summary(m2)##
## Call:
## lm(formula = phosphorus ~ arid.gradient, data = nutrient.mean)
##
## Residuals:
## 1 2 3 4 5 6 7
## 4.3035 1.1485 -2.0485 1.1930 -0.2791 -5.5020 1.1845
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.191 4.527 3.135 0.0258 *
## arid.gradient 1.780 1.298 1.371 0.2287
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.383 on 5 degrees of freedom
## Multiple R-squared: 0.2732, Adjusted R-squared: 0.1278
## F-statistic: 1.879 on 1 and 5 DF, p-value: 0.2287m3 <- lm(potassium~poly(arid.gradient,2), data=nutrient.mean)
summary(m3)##
## Call:
## lm(formula = potassium ~ poly(arid.gradient, 2), data = nutrient.mean)
##
## Residuals:
## 1 2 3 4 5 6 7
## -3.551 -65.927 32.611 57.492 -25.220 52.524 -47.931
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 186.56 22.69 8.221 0.00119 **
## poly(arid.gradient, 2)1 189.58 60.04 3.157 0.03427 *
## poly(arid.gradient, 2)2 190.31 60.04 3.170 0.03387 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 60.04 on 4 degrees of freedom
## Multiple R-squared: 0.8334, Adjusted R-squared: 0.7502
## F-statistic: 10.01 on 2 and 4 DF, p-value: 0.02774annual plant community repsonse
spp.data <- read.csv("Data/ERG.communitydata.csv")
spp.data[is.na(spp.data)] <- 0
mean.spp <- spp.data %>% group_by(Year,Site, Microsite) %>% summarize(abd=mean(Abundance), richness=mean(Richness), biomass=mean(Biomass))
mean.spp <- data.frame(mean.spp)
## collect aridity gradient data
site.climate <- rbind(s1.mean,s2.mean)
site.climate[,"Year"] <- c(rep("2016",7),rep("2017",7))
## combine climate data with community data
mean.spp <- merge(mean.spp, site.climate, by.y=c("Year","Site"))
# ## community responses season1
# mean.spp2016 <- subset(mean.spp, Year==2016)
#
# ## richness
# plot(mean.spp2016[mean.spp2016$Microsite=="open","arid.gradient"], mean.spp2016[mean.spp2016$Microsite=="open","richness"], ylab="average site richness", xlab="aridity gradient")
# points(mean.spp2016[mean.spp2016$Microsite=="shrub","arid.gradient"],mean.spp2016[mean.spp2016$Microsite=="shrub","richness"], pch=19)
#
# m1 <- lm(richness~arid.gradient*Microsite, data=mean.spp2016)
# summary(m1)
#
# ## abundance
# plot(mean.spp2016[mean.spp2016$Microsite=="open","arid.gradient"], mean.spp2016[mean.spp2016$Microsite=="open","abd"], ylab="average site abundance", xlab="aridity gradient")
# points(mean.spp2016[mean.spp2016$Microsite=="shrub","arid.gradient"],mean.spp2016[mean.spp2016$Microsite=="shrub","abd"], pch=19)
#
# m2 <- lm(abd~arid.gradient*Microsite, data=mean.spp2016)
# summary(m2)
#
# ## biomass
# plot(mean.spp2016[mean.spp2016$Microsite=="open","arid.gradient"], mean.spp2016[mean.spp2016$Microsite=="open","biomass"], ylab="average site biomass", xlab="aridity gradient")
# points(mean.spp2016[mean.spp2016$Microsite=="shrub","arid.gradient"],mean.spp2016[mean.spp2016$Microsite=="shrub","biomass"], pch=19)
#
# m3 <- lm(biomass~arid.gradient*Microsite, data=mean.spp2016)
# summary(m3)
#
# ##season2
#
# ## community responses season2
#
# mean.spp2017 <- subset(mean.spp, Year==2017)
#
#
# ## richness
# plot(mean.spp2017[mean.spp2017$Microsite=="open","arid.gradient"], mean.spp2017[mean.spp2017$Microsite=="open","richness"])
# points(mean.spp2017[mean.spp2017$Microsite=="shrub","arid.gradient"],mean.spp2017[mean.spp2017$Microsite=="shrub","richness"], pch=19)
#
# m1 <- lm(richness~arid.gradient*Microsite, data=mean.spp2017)
# summary(m1)
#
# ## abundance
# plot(mean.spp2017[mean.spp2017$Microsite=="open","arid.gradient"], mean.spp2017[mean.spp2017$Microsite=="open","abd"])
# points(mean.spp2017[mean.spp2017$Microsite=="shrub","arid.gradient"],mean.spp2017[mean.spp2017$Microsite=="shrub","abd"], pch=19)
#
# m2 <- lm(abd~arid.gradient*Microsite, data=mean.spp2017)
# summary(m2)
#
# ## biomass
# plot(mean.spp2017[mean.spp2017$Microsite=="open","arid.gradient"], mean.spp2017[mean.spp2017$Microsite=="open","biomass"])
# points(mean.spp2017[mean.spp2017$Microsite=="shrub","arid.gradient"],mean.spp2017[mean.spp2017$Microsite=="shrub","biomass"], pch=19)
#
# m3 <- lm(biomass~arid.gradient*Microsite, data=mean.spp2017)
# summary(m3)
## both seasons
## richness
plot(mean.spp[mean.spp$Microsite=="open","arid.gradient"], mean.spp[mean.spp$Microsite=="open","richness"], ylab="average site richness", xlab="aridity gradient")
points(mean.spp[mean.spp$Microsite=="shrub","arid.gradient"],mean.spp[mean.spp$Microsite=="shrub","richness"], pch=19) 
m1 <- lm(richness~arid.gradient*Microsite, data=mean.spp)
summary(m1)##
## Call:
## lm(formula = richness ~ arid.gradient * Microsite, data = mean.spp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.56169 -0.66225 -0.07639 0.66154 1.76433
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.9206 0.8381 5.871 4.68e-06 ***
## arid.gradient 0.8699 0.2730 3.186 0.00397 **
## Micrositeshrub -1.8968 1.1852 -1.600 0.12260
## arid.gradient:Micrositeshrub -0.5531 0.3861 -1.432 0.16490
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9061 on 24 degrees of freedom
## Multiple R-squared: 0.3356, Adjusted R-squared: 0.2526
## F-statistic: 4.042 on 3 and 24 DF, p-value: 0.01849## abundance
plot(mean.spp[mean.spp$Microsite=="open","arid.gradient"], mean.spp[mean.spp$Microsite=="open","abd"], ylab="average site abundance", xlab="aridity gradient")
points(mean.spp[mean.spp$Microsite=="shrub","arid.gradient"],mean.spp[mean.spp$Microsite=="shrub","abd"], pch=19) 
m2 <- lm(abd~arid.gradient*Microsite, data=mean.spp)
summary(m2)##
## Call:
## lm(formula = abd ~ arid.gradient * Microsite, data = mean.spp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -64.997 -26.611 -5.081 25.895 108.541
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 162.65 38.28 4.249 0.00028 ***
## arid.gradient 41.28 12.47 3.310 0.00294 **
## Micrositeshrub 39.22 54.14 0.724 0.47582
## arid.gradient:Micrositeshrub 10.04 17.64 0.569 0.57457
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 41.39 on 24 degrees of freedom
## Multiple R-squared: 0.5409, Adjusted R-squared: 0.4835
## F-statistic: 9.426 on 3 and 24 DF, p-value: 0.000268## biomass
plot(mean.spp[mean.spp$Microsite=="open","arid.gradient"], ihs(mean.spp[mean.spp$Microsite=="open","biomass"]), ylab="average site biomass", xlab="aridity gradient")
points(mean.spp[mean.spp$Microsite=="shrub","arid.gradient"],ihs(mean.spp[mean.spp$Microsite=="shrub","biomass"]), pch=19)
m3 <- lm(ihs(biomass)~arid.gradient*Microsite, data=mean.spp)
summary(m3)##
## Call:
## lm(formula = ihs(biomass) ~ arid.gradient * Microsite, data = mean.spp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6596 -0.5085 -0.0684 0.3314 1.2077
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.6251 0.5577 4.707 8.75e-05 ***
## arid.gradient 0.5873 0.1817 3.233 0.00355 **
## Micrositeshrub 1.6603 0.7888 2.105 0.04596 *
## arid.gradient:Micrositeshrub 0.3698 0.2569 1.439 0.16303
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.603 on 24 degrees of freedom
## Multiple R-squared: 0.6498, Adjusted R-squared: 0.606
## F-statistic: 14.84 on 3 and 24 DF, p-value: 1.128e-05cleaned up community response to gradient
RDA community and environmental data
community <- read.csv("Data/ERG.communitydata.csv")
community[is.na(community)] <- 0
##2016
## add enviro data to dataframe
comm2016 <- subset(community, Year==2016)
site.vars2016[,"Site"] <- as.factor(row.names(site.vars2016))
comm2016 <- merge(comm2016, site.vars2016, by="Site")
## transform spp data
community.trans <- decostand(comm2016[,13:53], "hellinger")
rda1 <- rda(community.trans, comm2016[,54:57])
(R2adj <- RsquareAdj(rda1)$adj.r.squared) ## 50.4% of variation explained## [1] 0.3871872## build byplot manually
par(mar=c(4.5,4.5,0.5,0.5))
plot(rda1, type="n", xlim=c(-1,1), ylim=c(-1.5,1.5))
## calculate priority
spp.priority <- colSums(comm2016[,13:53])
## plot RDA1
colvec <- rep(c("blue","dodgerblue3","cyan","yellow2","orange","tomato","red2"),each=60)
points(rda1, display = "sites", col = "black", pch = c(21), bg = colvec)
orditorp(rda1, display = "species", cex = 0.7, col = "darkred", priority=spp.priority, air=0.5)
text(rda1, display = "bp", col = "blue", cex = 0.8)
legend("bottomright", pch=c(21), legend=unique(comm2016$Site), pt.bg=unique(colvec), cex=1)
## ordikplot to adjust species locations
## 2017
## add enviro data to dataframe
comm2017 <- subset(community, Year==2017)
site.vars[,"Site"] <- as.factor(row.names(site.vars))
comm2017 <- merge(comm2017, site.vars, by="Site")
#
# ## transform spp data
# community.trans <- decostand(comm2017[,13:53], "hellinger")
#
# rda2 <- rda(community.trans, comm2017[,54:58])
# (R2adj <- RsquareAdj(rda2)$adj.r.squared) ## 50.4% of variation explained
#
#
# ## build byplot manually
# par(mar=c(4.5,4.5,0.5,0.5))
# plot(rda2, type="n", xlim=c(-1,1), ylim=c(-1.5,1.5))
#
# ## calculate priority
# spp.priority <- colSums(comm2017[,13:53])
#
# ## plot RDA2
# colvec <- rep(c("blue","dodgerblue3","cyan","yellow2","orange","tomato","red2"),each=60)
# points(rda2, display = "sites", col = "black", pch = c(21), bg = colvec)
# orditorp(rda2, display = "species", cex = 0.7, col = "darkred", priority=spp.priority, air=0.5)
# text(rda2, display = "bp", col = "blue", cex = 0.8)
# legend("bottomright", pch=c(21), legend=unique(comm2016$Site), pt.bg=unique(colvec), cex=1)
# ## ordikplot to adjust species locationsPhytometer analysis
census <- read.csv("Data/ERG.phytometer.census.csv")
census[is.na(census)] <- 0
census[,"phyto.abd"] <- rowSums(census[,c("Phacelia","Plantago","Salvia")])
## calculate the number of seeds per gram
seed.mass <- read.csv("Data/seed.mass.csv")
seed.mass.avg <- seed.mass %>% group_by(species) %>% summarize(avg.seed.gram=mean(seed.number),se.seed.gram=se(seed.number))
seed.mass.avg <- data.frame(seed.mass.avg)
seed.mass.avg## species avg.seed.gram se.seed.gram
## 1 Amsinkia 235.2 7.611833
## 2 Caulanthus 3019.0 8.916277
## 3 Lipidium 752.4 4.261455
## 4 Monolopia 737.8 7.831986
## 5 Phacelia 772.0 5.300943
## 6 Plantago 558.2 6.865858
## 7 Salvia 980.0 9.612492## get proportion of seed germinated per 0.3 grams of seed
census[,"Phacelia.prop"] <- census[,"Phacelia"]/(seed.mass.avg$avg.seed.gram[which(seed.mass.avg$species=="Phacelia")]*0.3)
census[,"Plantago.prop"] <- census[,"Plantago"]/(seed.mass.avg$avg.seed.gram[which(seed.mass.avg$species=="Plantago")]*0.3)
census[,"Salvia.prop"] <- census[,"Salvia"]/(seed.mass.avg$avg.seed.gram[which(seed.mass.avg$species=="Salvia")]*0.3)
## beginning
census.int <- subset(census, Census=="emergence")
## all species
m1 <- glm.nb(Phacelia~ Microsite * Nutrient * Site + Year, data=census.int)
anova(m1, test="Chisq") ## Site * Micro and Year significant## Analysis of Deviance Table
##
## Model: Negative Binomial(0.8108), link: log
##
## Response: Phacelia
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 839 1079.71
## Microsite 1 10.304 838 1069.41 0.001328 **
## Nutrient 1 1.234 837 1068.17 0.266538
## Site 6 159.651 831 908.52 < 2.2e-16 ***
## Year 1 31.424 830 877.10 2.073e-08 ***
## Microsite:Nutrient 1 0.046 829 877.05 0.829874
## Microsite:Site 6 21.032 823 856.02 0.001811 **
## Nutrient:Site 6 5.542 817 850.48 0.476389
## Microsite:Nutrient:Site 6 8.030 811 842.45 0.235879
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Phacelia
m1 <- glm.nb(Phacelia~ Microsite * Nutrient * Site + Year, data=census.int)
anova(m1, test="Chisq") ## Site * Micro and Year significant## Analysis of Deviance Table
##
## Model: Negative Binomial(0.8108), link: log
##
## Response: Phacelia
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 839 1079.71
## Microsite 1 10.304 838 1069.41 0.001328 **
## Nutrient 1 1.234 837 1068.17 0.266538
## Site 6 159.651 831 908.52 < 2.2e-16 ***
## Year 1 31.424 830 877.10 2.073e-08 ***
## Microsite:Nutrient 1 0.046 829 877.05 0.829874
## Microsite:Site 6 21.032 823 856.02 0.001811 **
## Nutrient:Site 6 5.542 817 850.48 0.476389
## Microsite:Nutrient:Site 6 8.030 811 842.45 0.235879
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Plantago
m2 <- glm.nb(Plantago~ Microsite * Nutrient * Site + Year, data=census.int)
anova(m2, test="Chisq") ## Site * Micro and Year significant## Analysis of Deviance Table
##
## Model: Negative Binomial(0.5335), link: log
##
## Response: Plantago
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 839 984.30
## Microsite 1 6.256 838 978.04 0.01238 *
## Nutrient 1 1.060 837 976.98 0.30332
## Site 6 162.687 831 814.30 < 2.2e-16 ***
## Year 1 4.138 830 810.16 0.04194 *
## Microsite:Nutrient 1 2.388 829 807.77 0.12224
## Microsite:Site 6 43.000 823 764.77 1.167e-07 ***
## Nutrient:Site 6 6.391 817 758.38 0.38088
## Microsite:Nutrient:Site 6 4.615 811 753.77 0.59403
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Salvia
m3 <- glm.nb(Salvia~ Microsite * Nutrient * Site + Year, data=census.int)
anova(m3, test="Chisq") ## Site * Micro and Year significant + (micro x nutrient)## Analysis of Deviance Table
##
## Model: Negative Binomial(0.8273), link: log
##
## Response: Salvia
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 839 1210.25
## Microsite 1 5.924 838 1204.33 0.0149386 *
## Nutrient 1 1.513 837 1202.82 0.2186678
## Site 6 211.625 831 991.19 < 2.2e-16 ***
## Year 1 17.463 830 973.73 2.929e-05 ***
## Microsite:Nutrient 1 7.824 829 965.90 0.0051548 **
## Microsite:Site 6 26.421 823 939.48 0.0001858 ***
## Nutrient:Site 6 4.587 817 934.90 0.5977624
## Microsite:Nutrient:Site 6 1.566 811 933.33 0.9549811
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1census.int.plot <- gather(census.int, species, abundance, Phacelia:Salvia)
##calculate confidence interval
census.plot<- census.int.plot %>% group_by(Microsite, species, Site, Year) %>% summarize(avg=mean(abundance),ci=se(abundance)*1.96)
ggplot(census.plot, aes(x=Microsite, y=avg, fill=species))+
geom_bar(position=position_dodge(), stat="identity")+
geom_errorbar(aes(ymin=avg-ci, ymax=avg+ci),
width=.2, # Width of the error bars
position=position_dodge(.9))+ scale_fill_brewer() +ylab("Abundance")+ theme_Publication() + facet_grid(Year~Site) 
## end season
census.end <- subset(census, Census=="end")
census.end[,"Year"] <- as.factor(census.end$Year)
## run full model for phyto biomass
m1 <- aov(log(phyto.biomass)~ Site * Microsite * Nutrient + Year, data=subset(census.end, phyto.biomass>0))
summary(m1)## Df Sum Sq Mean Sq F value Pr(>F)
## Site 6 330.8 55.13 32.164 < 2e-16 ***
## Microsite 1 10.9 10.88 6.347 0.012127 *
## Nutrient 1 30.3 30.25 17.651 3.24e-05 ***
## Year 1 161.2 161.21 94.051 < 2e-16 ***
## Site:Microsite 6 23.9 3.99 2.327 0.031940 *
## Site:Nutrient 6 41.9 6.98 4.073 0.000551 ***
## Microsite:Nutrient 1 1.7 1.73 1.011 0.315196
## Site:Microsite:Nutrient 6 7.9 1.32 0.768 0.595068
## Residuals 424 726.7 1.71
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(m1, "Microsite")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = log(phyto.biomass) ~ Site * Microsite * Nutrient + Year, data = subset(census.end, phyto.biomass > 0))
##
## $Microsite
## diff lwr upr p adj
## shrub-open 0.3096905 0.06750868 0.5518722 0.0123243## all species end of season
m2 <- glm.nb(phyto.abd~ Microsite * Nutrient * Site + Year, data=census.end)
anova(m2, test="Chisq") ## Site and microsite*nutrient significant## Analysis of Deviance Table
##
## Model: Negative Binomial(0.455), link: log
##
## Response: phyto.abd
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 839 1009.51
## Microsite 1 1.059 838 1008.45 0.30341
## Nutrient 1 2.527 837 1005.92 0.11194
## Site 6 135.420 831 870.50 < 2e-16 ***
## Year 1 3.412 830 867.09 0.06474 .
## Microsite:Nutrient 1 4.586 829 862.50 0.03224 *
## Microsite:Site 6 5.340 823 857.16 0.50104
## Nutrient:Site 6 5.169 817 852.00 0.52238
## Microsite:Nutrient:Site 6 5.683 811 846.31 0.45956
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## run full model for phacelia biomass
m3 <- aov(log(Phacelia.biomass)~ Site * Microsite * Nutrient + Year, data=subset(census.end, Phacelia.biomass>0))
summary(m3) ## site and microsite## Df Sum Sq Mean Sq F value Pr(>F)
## Site 6 89.7 14.95 8.274 4.95e-08 ***
## Microsite 1 21.2 21.20 11.738 0.00074 ***
## Nutrient 1 10.2 10.19 5.642 0.01845 *
## Year 1 95.5 95.49 52.857 7.42e-12 ***
## Site:Microsite 6 23.6 3.93 2.173 0.04692 *
## Site:Nutrient 6 8.0 1.34 0.742 0.61612
## Microsite:Nutrient 1 0.0 0.00 0.001 0.97552
## Site:Microsite:Nutrient 5 4.1 0.82 0.452 0.81122
## Residuals 205 370.3 1.81
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(m3, "Microsite")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = log(Phacelia.biomass) ~ Site * Microsite * Nutrient + Year, data = subset(census.end, Phacelia.biomass > 0))
##
## $Microsite
## diff lwr upr p adj
## shrub-open 0.6312547 0.2667053 0.9958041 0.0007714## run full model for plantago biomass
m4 <- aov(log(Plantago.biomass)~ Site * Microsite * Nutrient + Year, data=subset(census.end, Plantago.biomass>0))
summary(m4) ## site and year + sitexmicrosite## Df Sum Sq Mean Sq F value Pr(>F)
## Site 6 116.35 19.39 19.843 2.32e-16 ***
## Microsite 1 1.96 1.96 2.001 0.160
## Nutrient 1 32.03 32.03 32.778 6.80e-08 ***
## Year 1 138.38 138.38 141.597 < 2e-16 ***
## Site:Microsite 6 4.87 0.81 0.830 0.549
## Site:Nutrient 4 1.99 0.50 0.508 0.730
## Microsite:Nutrient 1 3.17 3.17 3.244 0.074 .
## Site:Microsite:Nutrient 3 0.65 0.22 0.220 0.882
## Residuals 130 127.05 0.98
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(m4, "Microsite")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = log(Plantago.biomass) ~ Site * Microsite * Nutrient + Year, data = subset(census.end, Plantago.biomass > 0))
##
## $Microsite
## diff lwr upr p adj
## shrub-open -0.2179629 -0.5348829 0.0989571 0.1759823## run full model for plantago biomass
m5 <- aov(log(Salvia.biomass)~ Site * Microsite * Nutrient + Year, data=subset(census.end, Salvia.biomass>0))
summary(m5) ## site and year + sitexmicrosite## Df Sum Sq Mean Sq F value Pr(>F)
## Site 6 170.1 28.35 19.982 < 2e-16 ***
## Microsite 1 0.0 0.01 0.008 0.929221
## Nutrient 1 21.4 21.43 15.103 0.000124 ***
## Year 1 44.1 44.05 31.045 5.39e-08 ***
## Site:Microsite 6 6.8 1.14 0.802 0.568988
## Site:Nutrient 6 46.2 7.70 5.424 2.35e-05 ***
## Microsite:Nutrient 1 0.5 0.55 0.387 0.534459
## Site:Microsite:Nutrient 6 12.9 2.16 1.520 0.171069
## Residuals 317 449.8 1.42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(m5, "Microsite")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = log(Salvia.biomass) ~ Site * Microsite * Nutrient + Year, data = subset(census.end, Salvia.biomass > 0))
##
## $Microsite
## diff lwr upr p adj
## shrub-open 0.01133515 -0.2406604 0.2633307 0.9295351## run full model for phacelia abundance
m6 <- glm.nb(Phacelia.biomass~ Site * Microsite * Nutrient + Year, data=subset(census.end, Phacelia>0))
summary(m3) ## site and microsite## Df Sum Sq Mean Sq F value Pr(>F)
## Site 6 89.7 14.95 8.274 4.95e-08 ***
## Microsite 1 21.2 21.20 11.738 0.00074 ***
## Nutrient 1 10.2 10.19 5.642 0.01845 *
## Year 1 95.5 95.49 52.857 7.42e-12 ***
## Site:Microsite 6 23.6 3.93 2.173 0.04692 *
## Site:Nutrient 6 8.0 1.34 0.742 0.61612
## Microsite:Nutrient 1 0.0 0.00 0.001 0.97552
## Site:Microsite:Nutrient 5 4.1 0.82 0.452 0.81122
## Residuals 205 370.3 1.81
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(m3, "Microsite")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = log(Phacelia.biomass) ~ Site * Microsite * Nutrient + Year, data = subset(census.end, Phacelia.biomass > 0))
##
## $Microsite
## diff lwr upr p adj
## shrub-open 0.6312547 0.2667053 0.9958041 0.0007714## run full model for plantago biomass
m7 <- aov(log(Plantago)~ Site * Microsite * Nutrient + Year, data=subset(census.end, Plantago>0))
summary(m4) ## site and year + sitexmicrosite## Df Sum Sq Mean Sq F value Pr(>F)
## Site 6 116.35 19.39 19.843 2.32e-16 ***
## Microsite 1 1.96 1.96 2.001 0.160
## Nutrient 1 32.03 32.03 32.778 6.80e-08 ***
## Year 1 138.38 138.38 141.597 < 2e-16 ***
## Site:Microsite 6 4.87 0.81 0.830 0.549
## Site:Nutrient 4 1.99 0.50 0.508 0.730
## Microsite:Nutrient 1 3.17 3.17 3.244 0.074 .
## Site:Microsite:Nutrient 3 0.65 0.22 0.220 0.882
## Residuals 130 127.05 0.98
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(m4, "Microsite")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = log(Plantago.biomass) ~ Site * Microsite * Nutrient + Year, data = subset(census.end, Plantago.biomass > 0))
##
## $Microsite
## diff lwr upr p adj
## shrub-open -0.2179629 -0.5348829 0.0989571 0.1759823## run full model for plantago biomass
m8 <- aov(log(Salvia)~ Site * Microsite * Nutrient + Year, data=subset(census.end, Salvia>0))
summary(m5) ## site and year + sitexmicrosite## Df Sum Sq Mean Sq F value Pr(>F)
## Site 6 170.1 28.35 19.982 < 2e-16 ***
## Microsite 1 0.0 0.01 0.008 0.929221
## Nutrient 1 21.4 21.43 15.103 0.000124 ***
## Year 1 44.1 44.05 31.045 5.39e-08 ***
## Site:Microsite 6 6.8 1.14 0.802 0.568988
## Site:Nutrient 6 46.2 7.70 5.424 2.35e-05 ***
## Microsite:Nutrient 1 0.5 0.55 0.387 0.534459
## Site:Microsite:Nutrient 6 12.9 2.16 1.520 0.171069
## Residuals 317 449.8 1.42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(m5, "Microsite")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = log(Salvia.biomass) ~ Site * Microsite * Nutrient + Year, data = subset(census.end, Salvia.biomass > 0))
##
## $Microsite
## diff lwr upr p adj
## shrub-open 0.01133515 -0.2406604 0.2633307 0.9295351census.end.plot <- gather(census.end, species, biomass, Phacelia.biomass:Salvia.biomass)
##calculate confidence interval
census.plot<- census.end.plot %>% group_by(Microsite, species, Site, Year) %>% summarize(avg=mean(biomass),ci=se(biomass)*1.96)
ggplot(census.plot, aes(x=Microsite, y=avg, fill=species))+
geom_bar(position=position_dodge(), stat="identity")+
geom_errorbar(aes(ymin=avg-ci, ymax=avg+ci),
width=.2, # Width of the error bars
position=position_dodge(.9))+ scale_fill_brewer() +ylab("Abundance")+ theme_Publication() + facet_grid(Year~Site) 
### RDA
#
# ##collect environmental variables for site
# nutrients <- read.csv("Data/ERG.soilnutrients.csv")
# nutrients.mean <- nutrients %>% group_by(gradient,site, microsite) %>% summarise_each(funs(mean))
# nutrients.vars <- data.frame(nutrients.mean)
#
# ## minimum explanation of variables
# ambient2016[is.na(ambient2016)] <- 0
# community <- ambient2016 %>% group_by(Gradient,Site, Microsite) %>% summarise_each(funs(sum))
# community <- data.frame(community)
#
# end.census <- subset(census2016, census=="end")
# end.census <- end.census %>% group_by(Gradient,Site, Microsite) %>% summarise_each(funs(mean))
# end.census<- data.frame(end.census)
#
# ## daily means
# day.mean <- aggregate(HOBO.data, by=list(day=HOBO.data$Day, micro=HOBO.data$Microsite, site=HOBO.data$Site), mean)
#
# ## summarize daily means across sites
# means <- day.mean %>% group_by(gradient, micro) %>% summarize(temp=mean(Temp),rh=mean(RH),temp.se=se(Temp),rh.se=se(RH))
# means <- data.frame(means)
#
#
# envs <- data.frame(swc=end.census[,"swc"],nutrients.vars[,c("N","P","K")], means[,c("temp","rh")])
#
# ## hellinger transformation
# community.trans <- decostand(community[,12:42], "hellinger")
#
# ## rda with environmental variables
# rda1 <- rda(community.trans, envs)
# anova(rda1)
# (R2adj <- RsquareAdj(rda1)$adj.r.squared) ## 25.7% of variation explained
#
#
# ## build byplot manually
# par(mar=c(4.5,4.5,0.5,0.5))
# plot(rda1, type="n", xlim=c(-1,1))
#
# with(community, levels(Site))
#
# spp.priority <- colSums(community[,12:42])
#
# colvec <- c("blue","dodgerblue3","cyan","yellow2","orange","tomato","red2")
# with(community, points(rda1, display = "sites", col = "black", pch = c(21,22), bg = colvec[Site]))
# orditorp(rda1, display = "species", cex = 0.7, col = "darkred", priority=spp.priority, air=1,)
# text(rda1, display = "bp", col = "blue", cex = 0.8)
# legend("bottomright", pch=c(21), legend=community$Site[community$Microsite=="shrub"], pt.bg=colvec, cex=1)
#
# ## check variable explanation
# vif.cca(rda1)
#
# v.clim <- cbind(envs[,c("temp","rh")])
# v.nut <- cbind(envs[,c("N","P","K")])
# v.swc <- envs[,"swc"]
#
# x <- varpart(community.trans,v.clim,v.nut,v.swc)
#
# showvarparts(3)Difference in effect size
library(bootES)## Loading required package: boot##
## Attaching package: 'boot'## The following object is masked from 'package:lattice':
##
## melanomaeffect.cal <- function(x){ bootES(x, R=999, data.col="Biomass", group.col="Microsite", contrast=c(shrub=1,open=-1), effect.type=c("cohens.d"))}
effect.site <- by(subset(comm2016, Site != "Barstow"), comm2016$Site[comm2016$Site != "Barstow"], FUN=effect.cal)
site.summary <- rbind(summary(effect.site$Panoche),summary(effect.site$Cuyama),summary(effect.site$TejonRanch),data.frame(stat=0,ci.low=0,ci.high=0,bias=0,std.error=0),summary(effect.site$MojavePreserve),summary(effect.site$SheepholeValley),summary(effect.site$Tecopa))
site.summary[5,] <- 0
site.summary <- sapply(site.summary, as.numeric)
par(mar=c(4.5,4.5,.5,.5))
plot(c(1:7),site.summary[,"stat"], xlim=c(0.5,7.5), ylim=c(-1,3), pch=19, cex=1.4, ylab="Hedge's G", xlab="", xaxt="n")
arrows(c(1:7),as.numeric(site.summary[,"ci.high"]),c(1:7),site.summary[,"ci.low"], angle=90, code=3, length=0, lwd=2)
abline(h=0, lty=2, lwd=2)
axis(1, c(1:7), labels=site.names, cex.axis=1)
## phytometer
effect.cal <- function(x){ bootES(x, R=999, data.col="phyto.biomass", group.col="Microsite", contrast=c(shrub=1,open=-1), effect.type=c("cohens.d"))}
census2016 <- subset(census, Census=="end" & Year ==2016)
effect.site <- by(subset(census2016, Site != "Barstow"), census2016$Site[census2016$Site != "Barstow"], FUN=effect.cal)
site.summary <- rbind(summary(effect.site$Panoche),summary(effect.site$Cuyama),summary(effect.site$TejonRanch),data.frame(stat=0,ci.low=0,ci.high=0,bias=0,std.error=0),summary(effect.site$MojavePreserve),summary(effect.site$SheepholeValley),summary(effect.site$Tecopa))
site.summary <- sapply(site.summary, as.numeric)
par(mar=c(4.5,4.5,.5,.5))
plot(c(1:7),site.summary[,"stat"], xlim=c(0.5,7.5), ylim=c(-1,3), pch=19, cex=1.4, ylab="Hedge's G", xlab="", xaxt="n")
arrows(c(1:7),site.summary[,"ci.high"],c(1:7),site.summary[,"ci.low"], angle=90, code=3, length=0, lwd=2)
abline(h=0, lty=2, lwd=2)
axis(1, c(1:7), labels=site.names, cex.axis=1)
Rii between years
rii.dat <- rii(community, j=1:8, var=c("Abundance","Richness","Biomass"))
rii.mean <- rii.dat %>% group_by(Year, Gradient, Site) %>% summarize(avg=mean(Biomass),se=se(Biomass))
rii.mean <- data.frame(rii.mean)
plot(rii.mean[rii.mean$Year==2016,"Gradient"]-.1,rii.mean[rii.mean$Year==2016,"avg"], ylim=c(-1,1), pch=19, cex=1.5, xlim=c(0.5,7.5), ylab="Rii Biomass", xlab="Site", xaxt="n", cex.lab=1.5)
axis(1, 1:7, lab=unique(rii.mean$Site))
error.bar(rii.mean[rii.mean$Year==2016,"Gradient"]-.1,rii.mean[rii.mean$Year==2016,"avg"],rii.mean[rii.mean$Year==2016,"se"])
error.bar(rii.mean[rii.mean$Year==2017,"Gradient"]+.1,rii.mean[rii.mean$Year==2017,"avg"],rii.mean[rii.mean$Year==2017,"se"])
points(rii.mean[rii.mean$Year==2017,"Gradient"]+.1,rii.mean[rii.mean$Year==2017,"avg"], pch=21, bg="Grey50", cex=1.5)
abline(h=0, lwd=2, lty=2)
legend(6.3,-0.55, c("2016","2017"), pch=22, pt.bg=c("Black","Grey50"), cex=1.6)
## add year column to precipitation data
precip[,"Year"] <- ifelse(precip$season=="season.1","2016","2017")
rii.precip <- merge(precip,rii.mean, by=c("Year","Gradient"))